
% function R = romberg(f, a, b, n)
% format long

a = 0.89;
b = 0.6;
n = 5;
f = sin(a);

for i = 1 : n
    h = (b - a) / 2^i;
    s = 0;
    for k = 1 : 2^(i-1),
        s = s + feval(f, a + (2*k - 1)*h);
    end
    R(i+1, 0+1) = R(i-1+1, 0+1)/2 + h*s;
end